# install.packages("Seurat")
# install.packages("dotCall64")
# Seurat v4 was just loaded with SeuratObject v5; disabling v5 assays and validation
# routines, and ensuring assays work in strict v3/v4 compatibility mode
# 这段文字描述了在加载Seurat v4版本时，同时加载了SeuratObject v5对象，并禁用了一些特定的实验方法和验证例程，以确保在严格遵循v3/v4兼容性模式下进行数据分析。在加载Seurat v4版本时，同时加载了SeuratObject v5对象，并禁用了一些特定的实验方法和验证例程，以确保在严格遵循v3/v4兼容性模式下进行数据分析。
library(Seurat)
library(ggplot2)
library(cowplot)
# install.packages("scater")
# install.packages("BiocManager")
# BiocManager::install("scater")
library(scater)
library(scran)
library(BiocParallel)
library(BiocNeighbors)
library(data.table)
library(dplyr)
library(Matrix)
# BiocManager::install("clustree")
library(clustree)
library(pheatmap)
setwd("C:/Users/forbing36/Desktop/单细胞生信/GSE212966")


#单组分析####
#**********************************************
# 需要把数据放到不同文件目录下，且文件名应该是
# barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz
# 没有gz后缀也可以
#**********************************************

##ADJ1
a1 <- Read10X(data.dir = "./data/input/GSE212966_RAW/ADJ1")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="ADJ1.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "ADJ1.rds")


##ADJ2
a1 <- Read10X(data.dir = "./GSE212966_RAW/ADJ2")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="ADJ2.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "ADJ2.rds")

##ADJ3
a1 <- Read10X(data.dir = "./GSE212966_RAW/ADJ3")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="ADJ3.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "ADJ3.rds")

##ADJ4
a1 <- Read10X(data.dir = "./GSE212966_RAW/ADJ4")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="ADJ4.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "ADJ4.rds")

##ADJ5
a1 <- Read10X(data.dir = "./GSE212966_RAW/ADJ5")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="ADJ5.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "ADJ5.rds")

##PDAC1
a1 <- Read10X(data.dir = "./GSE212966_RAW/PDAC1")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="PDAC1.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "PDAC1.rds")


##PDAC2
a1 <- Read10X(data.dir = "./GSE212966_RAW/PDAC2")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="PDAC2.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "PDAC2.rds")

##PDAC3
a1 <- Read10X(data.dir = "./GSE212966_RAW/PDAC3")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="PDAC3.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "PDAC3.rds")

##PDAC4
a1 <- Read10X(data.dir = "./GSE212966_RAW/PDAC4")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="PDAC4.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "PDAC4.rds")


##PDAC5
a1 <- Read10X(data.dir = "./GSE212966_RAW/PDAC5")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="PDAC5.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "PDAC5.rds")



##PDAC6
a1 <- Read10X(data.dir = "./GSE212966_RAW/PDAC6")
pbmc <- CreateSeuratObject(counts = a1, project = "a1", min.cells =3, min.features=200)
pbmc
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pdf(file="PDAC6.pdf", width = 13, height = 6)
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

#pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#Determine the ‘dimensionality’ of the dataset
#pbmc <- JackStraw(pbmc, num.replicate = 100,dims = 40)
#pbmc <- ScoreJackStraw(pbmc, dims = 1:40)
#JackStrawPlot(pbmc, dims = 1:40)
#ElbowPlot(pbmc,ndims = 40)
#確定下面的dims
pbmc <- RunUMAP(pbmc, dims = 1:20)
saveRDS(pbmc, file = "PDAC6.rds")


# 总体分析####
##1. 读取分组数据####
# 使用for循环创建变量并读入数据
num_vars <- 6
for (i in 1:num_vars) {
  if(i>=3 && i<=5){next}  #ADJ3，ADJ4，ADJ5数据缺失，跳过执行
  var_name <- paste("ADJ", i, sep = "")
  assign(var_name,readRDS(file=paste0(var_name,".rds")))# 为变量赋值
}
for (i in 1:num_vars) {
  var_name <- paste("PDAC", i, sep = "")
  assign(var_name,readRDS(file=paste0(var_name,".rds")))# 为变量赋值
}

##2. 标注分组信息####
ADJ1<-RenameCells(ADJ1,add.cell.id="ADJ1",for.merge=T)
ADJ1@meta.data$tech<-"ADJ"
ADJ1@meta.data$celltype<-"ADJ_ADJ1"

ADJ2<-RenameCells(ADJ2,add.cell.id="ADJ2",for.merge=T)
ADJ2@meta.data$tech<-"ADJ"
ADJ2@meta.data$celltype<-"ADJ_ADJ2"

ADJ3<-RenameCells(ADJ3,add.cell.id="ADJ3",for.merge=T)
ADJ3@meta.data$tech<-"ADJ"
ADJ3@meta.data$celltype<-"ADJ_ADJ3"

ADJ4<-RenameCells(ADJ4,add.cell.id="ADJ4",for.merge=T)
ADJ4@meta.data$tech<-"ADJ"
ADJ4@meta.data$celltype<-"ADJ_ADJ4"

ADJ5<-RenameCells(ADJ5,add.cell.id="ADJ5",for.merge=T)
ADJ5@meta.data$tech<-"ADJ"
ADJ5@meta.data$celltype<-"ADJ_ADJ5"

ADJ6<-RenameCells(ADJ6,add.cell.id="ADJ6",for.merge=T)
ADJ6@meta.data$tech<-"ADJ"
ADJ6@meta.data$celltype<-"ADJ_ADJ6"


PDAC1<-RenameCells(PDAC1,add.cell.id="PDAC1",for.merge=T)
PDAC1@meta.data$tech<-"PDAC"
PDAC1@meta.data$celltype<-"PDAC_PDAC1"

PDAC2<-RenameCells(PDAC2,add.cell.id="PDAC2",for.merge=T)
PDAC2@meta.data$tech<-"PDAC"
PDAC2@meta.data$celltype<-"PDAC_PDAC2"

PDAC3<-RenameCells(PDAC3,add.cell.id="PDAC3",for.merge=T)
PDAC3@meta.data$tech<-"PDAC"
PDAC3@meta.data$celltype<-"PDAC_PDAC3"

PDAC4<-RenameCells(PDAC4,add.cell.id="PDAC4",for.merge=T)
PDAC4@meta.data$tech<-"PDAC"
PDAC4@meta.data$celltype<-"PDAC_PDAC4"

PDAC5<-RenameCells(PDAC5,add.cell.id="PDAC5",for.merge=T)
PDAC5@meta.data$tech<-"PDAC"
PDAC5@meta.data$celltype<-"PDAC_PDAC5"

PDAC6<-RenameCells(PDAC6,add.cell.id="PDAC6",for.merge=T)
PDAC6@meta.data$tech<-"PDAC"
PDAC6@meta.data$celltype<-"PDAC_PDAC6"


##3. 数据合并####
##ADJ
# ADJ12<-merge(ADJ1,ADJ2)
# ADJ34<-merge(ADJ3,ADJ4)
# ADJ56<-merge(ADJ5,ADJ6)
# ADJ1234 <- merge(ADJ12,ADJ34)
# ADJ123456 <- merge(ADJ1234,ADJ56)

##PDAC
# PDAC12<-merge(PDAC1,PDAC2)
# PDAC34<-merge(PDAC3,PDAC4)
# PDAC56<-merge(PDAC5,PDAC6)
# PDAC1234 <- merge(PDAC12,PDAC34)
# PDAC123456 <- merge(PDAC1234,PDAC56)

####法2
# 创建一个包含代码的字符串
# 使用eval()函数执行代码
# 或者,将代码保存到一个文件中，然后使用source_file()函数执行
# writeLines(c('print("Hello, World!")'), "my_code.R")
# source_file("my_code.R")

#ADJ
code_str <- ""
for(i in 1:6){
  code_str <- paste(code_str,paste0(paste0("merge(ADJ,ADJ",i),")"),sep = " \n ")
}
code_str  # " \n merge(ADJ,ADJ1) \n merge(ADJ,ADJ2) \n merge(ADJ,ADJ3) \n merge(ADJ,ADJ4) \n merge(ADJ,ADJ5) \n merge(ADJ,ADJ6)"
eval(parse(text = code_str))

#PDAC🌟
code_str <- ""
for(i in 1:6){
  code_str <- paste(code_str,paste0(paste0("merge(PDAC,PDAC",i),")"),sep = " \n ")
}
code_str  # " \n merge(PDAC,PDAC1) \n merge(PDAC,PDAC2) \n merge(PDAC,PDAC3) \n merge(PDAC,PDAC4) \n merge(PDAC,PDAC5) \n merge(PDAC,PDAC6)"
eval(parse(text = code_str))



## 4. before integrate####
library(Seurat)
library(ggplot2)
hms<-readRDS(file="ADJ_PDAC_before_integrate.rds")

###ALL_Feature_Count_Percent_VlinPlot####
hms[["percent.mt"]] <- PercentageFeatureSet(hms, pattern = "^MT-")
dim(hms[["percent.mt"]])#全为0的话是不对，修改上面的pattern正则表达
# pdf(file="ADJ_PDAC_VlnPlot.pdf", width = 13, height = 6)
pdf(file="ALL_Feature_Count_Percent_VlinPlot.pdf", width = 13, height = 6)
VlnPlot(hms, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,pt.size = 0)
plot1 <- FeatureScatter(hms, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(hms, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
dev.off()

# DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
# DimHeatmap(pbmc, dims = 1:30, cells = 500, balanced = TRUE)


# ？？？
pancreas <- NormalizeData(object = hms, normalization.method = "LogNormalize", scale.factor = 1e4)
pancreas <- FindVariableFeatures(pancreas, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
pancreas <- ScaleData(pancreas, verbose = FALSE)
pancreas <- RunPCA(pancreas, npcs = 30, verbose = FALSE)
pancreas <- RunUMAP(pancreas, reduction = "pca", dims = 1:30)
p1 <- DimPlot(pancreas, reduction = "umap", group.by = "tech")
p2 <- DimPlot(pancreas, reduction = "umap", group.by = "celltype", label = TRUE, repel = TRUE) + 
  NoLegend()
###hms_before_integrate_DimPlot####
# pdf(file="ADJ_PDAC_DimPlot.pdf", width = 13, height = 6)
pdf(file="hms_before_integrate_DimPlot.pdf", width = 13, height = 6)
plot_grid(p1,p2)
dev.off()

##5. after integrate####
pancreas.list <- SplitObject(pancreas, split.by = "celltype")
for (i in 1: length(pancreas.list)) {
  pancreas.list[[i]] <- NormalizeData(pancreas.list[[i]], verbose = FALSE)
  pancreas.list[[i]] <- FindVariableFeatures(pancreas.list[[i]], selection.method = "vst", nfeatures = 2000, 
                                             verbose = FALSE)
}
#reference.list <- pancreas.list[c("Blood_P57","Tumor_P57","Blood_P58","Tumor_P58", 
#"Blood_P60","Tumor_P60", "Blood_P61", "Tumor_P61","Juxta_P60","Juxta_P61")]
reference.list <- pancreas.list[c("ADJ_ADJ1","ADJ_ADJ2","ADJ_ADJ3","ADJ_ADJ4","ADJ_ADJ5","ADJ_ADJ6",
                                  "PDAC_PDAC1","PDAC_PDAC2", "PDAC_PDAC3", "PDAC_PDAC4", "PDAC_PDAC5", "PDAC_PDAC6")]
pancreas.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)
pancreas.integrated <- IntegrateData(anchorset = pancreas.anchors, dims = 1:30)
DefaultAssay(pancreas.integrated) <- "integrated"
pancreas.integrated <- ScaleData(pancreas.integrated, verbose = FALSE)
pancreas.integrated <- RunPCA(pancreas.integrated, npcs = 30, verbose = FALSE)

print(pancreas.integrated[["pca"]], dims = 1:30, nfeatures = 5)
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(hms_cluster, reduction = "pca")

pancreas.integrated <- RunUMAP(pancreas.integrated, reduction = "pca", dims = 1:30)
###hms_after_integrated_DimPlot####
pdf(file="hms_after_integrated_DimPlot.pdf", width = 13, height = 6)
p1 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "tech")
p2 <- DimPlot(pancreas.integrated, reduction = "umap", group.by = "celltype")
plot_grid(p1,p2)
dev.off()
saveRDS(pancreas.integrated, file = "hms_after_integrated.rds")
# saveRDS(pancreas.integrated, file = "DAC_PDAC_after_integrated.rds")

#聚类获得亚群####
hms_individual_integrated<-readRDS(file="hms_after_integrated.rds")
# p1 <- DimPlot(hms_individual_integrated, reduction = "umap", group.by = "celltype")
# p1

##1. 确定dim参数(可以不出图)####
pbmc <- JackStraw(hms_individual_integrated, num.replicate = 100,dims = 40)
pbmc <- ScoreJackStraw(pbmc, dims = 1:30)
###hms_individual_integrated_JackStrawPlot_ElbowPlot####
pdf("hms_individual_integrated_JackStrawPlot_ElbowPlot.pdf",width = 13,height = 6)
JackStrawPlot(pbmc, dims = 1:30)
ElbowPlot(pbmc,ndims = 30)
# find how many 15cluster
ElbowPlot(hms_individual_integrated)
dev.off()

##2. 亚群聚类####
###2.1 亚组分群####
hms_neighbor<- FindNeighbors(hms_individual_integrated, dims = 1:30)
obj <- FindClusters(hms_neighbor, resolution = seq(0.5,1.2,by=0.1))
#resolution設置在0.4-1.2之間,越大clusters越多,查看拐點
####hms_neighbor_Clustree####
pdf("hms_neighbor_Clustree.pdf",width = 13,height = 6)
clustree(obj)
dev.off()
hms_cluster <- FindClusters( hms_neighbor, resolution = 1.2)
head(Idents(hms_cluster), 5)
hms_cluster<- RunUMAP(hms_cluster, dims = 1:30)
####hms_cluster_DimPlot####
pdf("hms_cluster_DimPlot.pdf",width = 13,height = 6)
DimPlot(hms_cluster, reduction = "umap")
dev.off()
saveRDS(hms_cluster, file = "hms_cluster_test_1.2.rds")
###2.2 确定markers####
scRNA.markers <- FindAllMarkers(hms_cluster, 
                                only.pos = TRUE,  #特异性高表达marker
                                min.pct = 0.05, 
                                logfc.threshold = 0.05
)
write.table(scRNA.markers,file="cellMarkers.txt",sep="\t",row.names=F,quote=F)

####top20_cell_markers####
#挑选每个细胞亚群中特意高表达的20个基因
top20 <- scRNA.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC) 
write.csv(file="top20_cell_markers.csv",top20)
DoHeatmap(hms_cluster, features = top20$gene) + NoLegend()

#整理成表格，只显示基因名字
top20_table=unstack(top20, gene ~ cluster)
names(top20_table)=gsub("X","cluster",names(top20_table))
write.csv(file="top20_marker_genes_1.2.csv",top20_table,row.names=F)


slotNames(hms_cluster) # 细胞及细胞中基因与RNA数量
hms_cluster@assays     # assay
dim(hms_cluster@meta.data)
View(hms_cluster@meta.data)



hms_cluster<-readRDS(file="./data/temp/hms_cluster_test_1.2.rds")
DimPlot(hms_cluster, reduction = "umap")
##3. 自动注释(⚠️不使用)####
# SingleR https://zhuanlan.zhihu.com/p/448269413
library(SingleR)
library(Seurat)
library(pheatmap)
# load("HumanPrimaryCellAtlas_hpca.se_human.RData")
load("BlueprintEncode_bpe.se_human.RData")

#进行SingleR注释
hms_cluster_for_SingleR <- GetAssayData(hms_cluster, slot="data") ##获取标准化矩阵
hms_cluster.hesc <- SingleR(test = hms_cluster_for_SingleR, ref = bpe.se, labels = bpe.se$label.main) #
hms_cluster.hesc
hms_cluster.hesc$labels

#seurat 和 SingleR的table表
a <- table(hms_cluster.hesc$labels, hms_cluster@meta.data$seurat_clusters)# pbmc@meta.data  pbmc的meta文件，包含了seurat的聚类结果
write.csv(a, "SingleR_by_pbe.csv")

#保存亚群名到hms_cluster
hms_cluster@meta.data$labels <-hms_cluster.hesc$labels
####sub_clusters_singleR_by_bpe####
pdf("hms_cluster_SingleR_by_bpe.pdf",width = 21,height = 8)
DimPlot(hms_cluster, group.by = c("seurat_clusters", "labels"),reduction = "umap",label = T)
dev.off()
saveRDS(hms_cluster, file = "hms_cluster_by_SingleR_labels.rds")

##3.手动注释####
hms_cluster <- readRDS(file="./data/temp/hms_cluster_by_SingleR_labels.rds")  

# markers1 <- c(
#   "COL1A1", "COL1A2", "COL3A1", "DCN",# Fibroblast
#   "CD79A","MS4A1", "LINC01781", "CD83",# B cell "CD19",
#   "IL7R", "CCR7", "LTB","GPR183",   # CD4+T-cell
#   "CLDN4", "KRT19", "KRT7", "MMP7", # Pancreatic ductal cell
#   "GZMA","CCL5", "GZMK", "NKG7",# CD8+T-cells
#   "CD68", "CD14", "CD163", # Macrophage
#   "VWF", "KDR", "PECAM1", "CLDN5"# Endothelial cell
# )
# markers2 <- c(
#   "PRSS1","CTRB1","PRSS2","REG1B",# Acinar cell
#   "ADIRF", "RGS5", "COL4A1", "PDGFRB",  # Pancreatic stellate cell
#   "KIT", "CPA3", "TPSAB1","TPSB2",  # Mast cell
#   "CHGA", "INS", "CHGB", "TTR", # Endocrine cell
#   "GNLY", "GZMB", "KLRD1", "PRF1",# NK cells
#   "IGHA1", "JCHAIN", "MZB1",# Plasma cell
#   "MKI67", "CDK1", "HMGB2", "STMN1",# Erythrocytes
#   "CDH19", "GPM6B"# Malignant cell
# )
markers1 <- c(
  "COL1A1", "COL1A2", "COL3A1", "DCN",# Fibroblast
  "CD79A","MS4A1", "LINC01781", "CD83",# B cell "CD19",
  "IL7R", "CCR7", "LTB","GPR183",   # T-cell 待定
  "CLDN4", "KRT19", "KRT7", "MMP7", # Pancreatic ductal cell
  "CD68", "CD14", "CD163", # Macrophage
  "VWF", "KDR", "PECAM1", "CLDN5"# Endothelial cell
)
markers2 <- c(
  "PRSS1","CTRB1","PRSS2","REG1B",# Acinar cell
  "ADIRF", "RGS5", "COL4A1", "PDGFRB",  # Pancreatic stellate cell
  "KIT", "CPA3", "TPSAB1","TPSB2",  # Mast cell
  "CHGA", "INS", "CHGB", "TTR", # Endocrine cell
  "GNLY", "GZMB", "KLRD1", "PRF1",# NK cells
  "IGHA1", "JCHAIN", "MZB1",# Plasma cell
  # "MKI67", "CDK1", "HMGB2", "STMN1",# Erythrocytes
  "CDH19", "LGI4","S100B","SCN7A")# Schwann cell
dp1 <-DotPlot(hms_cluster, features = markers1, dot.scale = 8) + RotatedAxis()
dp2 <-DotPlot(hms_cluster, features = markers2, dot.scale = 8) + RotatedAxis()
dp1 + dp2
####hms_cluster_手动标注####
pdf("./data/output/hms_cluster_标注验证.pdf",width = 20,height = 12)
dp1 + dp2
dev.off()

#test
markers3 <- c("CD3","CD3D","CD3E")
dp1 <- DotPlot(hms_cluster, features = markers3, dot.scale = 8) + RotatedAxis()
dp1
pdf("./data/output/1.pdf",width = 20,height = 12)
dp1
dev.off()
##4. 注释验证####
hms_cluster<-readRDS("./data/temp/hms_cluster_id_test.rds")
###4.1 方法1: DotPlot####
markers1 <- c(
  "COL1A1", "COL1A2", "COL3A1",# Fibroblast "CTHRC1",
  "MS4A1", "CD79A","LINC01781",# B cell "CD19", "CD83",
  "IL7R", "GZMK","TIGIT",   # T-cell 待定"BATF", 
  "CLDN4", "KRT19", "KRT7",  # Pancreatic ductal cell"MMP7",
  "CD68", "CD14", "CD163", # Macrophage
  "VWF",  "PECAM1", "CLDN5",# Endothelial cell"KDR",
  "GNLY", "GZMB", "KLRD1",# NK cells "PRF1",
  "PRSS1","CTRB1","PRSS2",# Acinar cell"REG1B",
  "G0S2", "S100A9", "S100A8",# Neitrophil "CXCL8",
  "ADIRF", "RGS5", "PDGFRB",  # Pancreatic stellate cell "COL4A1", 
  "CPA3", "TPSAB1","TPSB2",  # Mast cell"KIT", 
  "CHGA", "INS", "CHGB", # Endocrine cell"TTR", 
  "IGHA1", "JCHAIN", "MZB1",# Plasma cell
  # "MKI67", "CDK1", "HMGB2", "STMN1",# Erythrocytes
  "CDH19", "LGI4","S100B")# Schwann cell,"SCN7A"
dp1 <-DotPlot(hms_cluster, features = markers1, dot.scale = 5) + RotatedAxis()
dp1
#####DotPlot_全聚类_标注验证####
pdf("./data/output/DotPlot_全聚类_标注验证.pdf",width = 13,height = 6)
dp1
dev.off()


###4.2 方法2: FeaturePlot####
fp1 <- FeaturePlot(hms_cluster,reduction ="umap",features =c("HBB"))
fp1
###4.2 方法3: VlnPlot####
fp2 <- VlnPlot(hms_cluster, features = c("MS4A1", "CD79A"))
fp2

##5. 注释基因####
new.cluster.ids <- c("Fibroblast","B cell","T cell","Ductal cell","T cell",
                     "Macrophage","Fibroblast","Endothelial cell","Ductal cell",
                     "T cell","NK cell","Ductal cell","Ductal cell","Acinar cell",
                     "Neutrophil","Fibroblast","T cell","T cell","T cell","T cell",
                     "Stellate cell","Mast cell","T cell","Endocrine cell","T cell",
                     "Plasma cell","Stellate cell","Ductal cell","Plasma cell",
                     "T cell","Schwann cell","T cell","Endothelial cell",
                     "Endocrine cell","Macrophage","Fibroblast","Fibroblast","B cell")

names(new.cluster.ids) <- levels(hms_cluster)
#在hms_cluster@active.ident添加分类标记
hms_cluster_id<- RenameIdents(hms_cluster, new.cluster.ids)
hms_cluster_id@active.ident

####hms_cluster_id_DimPlot####
pdf("./data/output/hms_cluster_id_DimPlot1.pdf",width =23 ,height = 10)
dp1 <- DimPlot(hms_cluster_id, reduction = "umap", label = TRUE, pt.size = 0.5) 
dp2 <- DimPlot(hms_cluster_id, reduction = "umap", label = FALSE, pt.size = 0.5)
dp1+dp2
dev.off()
saveRDS(hms_cluster_id, file = "./data/temp/hms_cluster_id_test.rds")


##6. 分群绘图####
hms_cluster_id<-readRDS("./data/temp/hms_cluster_id_test.rds")

Fibroblast<-subset(hms_cluster_id, idents=c('Fibroblast'))
pdf("./data/output/Fibroblast_DimPlot.pdf",width = 8,height = 6)
DimPlot(Fibroblast, reduction = "umap")
dev.off()
saveRDS(Fibroblast, file="./data/temp/Fibroblast.rds")

B<-subset(hms_cluster_id, idents=c('B cell'))
pdf("./data/output/B_DimPlot.pdf",width = 8,height = 6)
DimPlot(B, reduction = "umap")
dev.off()
saveRDS(B, file="./data/temp/B.rds")

T<-subset(hms_cluster_id, idents=c('T cell'))
###T_DimPlot####
pdf("./data/output/T_DimPlot.pdf",width = 8,height = 6)
DimPlot(T, reduction = "umap")
dev.off()
saveRDS(T, file="./data/temp/T.rds")

Ductal<-subset(hms_cluster_id, idents=c('Ductal cell'))
pdf("./data/output/Ductal_DimPlot.pdf",width = 8,height = 6)
DimPlot(Ductal, reduction = "umap")
dev.off()
saveRDS(Ductal, file="./data/temp/Ductal.rds")

Macrophage<-subset(hms_cluster_id, idents=c('Macrophage'))
pdf("./data/output/Macrophage_DimPlot.pdf",width = 8,height = 6)
DimPlot(Macrophage, reduction = "umap")
dev.off()
saveRDS(Macrophage, file="./data/temp/Macrophage.rds")

Endothelial<-subset(hms_cluster_id, idents=c('Endothelial cell'))
pdf("./data/output/Endothelial_DimPlot.pdf",width = 8,height = 6)
DimPlot(Endothelial, reduction = "umap")
dev.off()
saveRDS(Endothelial, file="./data/temp/Endothelial.rds")

NK<-subset(hms_cluster_id, idents=c('NK cell'))
pdf("./data/output/NK_DimPlot.pdf",width = 8,height = 6)
DimPlot(NK, reduction = "umap")
dev.off()
saveRDS(NK, file="./data/temp/NK.rds")

Acinar<-subset(hms_cluster_id, idents=c('Acinar cell'))
pdf("./data/output/Acinar_DimPlot.pdf",width = 8,height = 6)
DimPlot(Acinar, reduction = "umap")
dev.off()
saveRDS(Acinar, file="./data/temp/Acinar.rds")

Neutrophil<-subset(hms_cluster_id, idents=c('Neutrophil'))
pdf("./data/output/Neutrophil_DimPlot.pdf",width = 8,height = 6)
DimPlot(Neutrophil, reduction = "umap")
dev.off()
saveRDS(Neutrophil, file="./data/temp/NK.rds")

Stellate<-subset(hms_cluster_id, idents=c('Stellate cell'))
pdf("./data/output/Stellate_DimPlot.pdf",width = 8,height = 6)
DimPlot(Stellate, reduction = "umap")
dev.off()
saveRDS(Stellate, file="./data/temp/Stellate.rds")

Mast<-subset(hms_cluster_id, idents=c('Mast cell'))
pdf("./data/output/Mast_DimPlot.pdf",width = 8,height = 6)
DimPlot(Mast, reduction = "umap")
dev.off()
saveRDS(Mast, file="./data/temp/Mast.rds")

Endocrine<-subset(hms_cluster_id, idents=c('Endocrine cell'))
pdf("./data/output/Endocrine_DimPlot.pdf",width = 8,height = 6)
DimPlot(Endocrine, reduction = "umap")
dev.off()
saveRDS(Endocrine, file="./data/temp/Endocrine.rds")

Plasma<-subset(hms_cluster_id, idents=c('Plasma cell'))
pdf("./data/output/Plasma_DimPlot.pdf",width = 8,height = 6)
DimPlot(Plasma, reduction = "umap")
dev.off()
saveRDS(Plasma, file="./data/temp/Plasma.rds")

Schwann<-subset(hms_cluster_id, idents=c('Schwann cell'))
pdf("./data/output/Schwann_DimPlot.pdf",width = 8,height = 6)
DimPlot(Schwann, reduction = "umap")
dev.off()
saveRDS(Schwann, file="./data/temp/Schwann.rds")


###展示亚群对比####
##input each cluster
Fibroblast <- readRDS(file = "./data/temp/Fibroblast.rds")
B <- readRDS(file = "./data/temp/B.rds")
T <- readRDS(file = "./data/temp/T.rds")
Ductal <- readRDS(file = "./data/temp/Ductal.rds")
Macrophage <- readRDS(file = "./data/temp/Macrophage.rds")
Endothelial <- readRDS(file = "./data/temp/Endothelial.rds")
NK <- readRDS(file = "./data/temp/NK.rds")
Acinar <- readRDS(file = "./data/temp/Acinar.rds")
Neutrophil <- readRDS(file = "./data/temp/NK.rds")
Stellate <- readRDS(file = "./data/temp/Stellate.rds")
Mast <- readRDS(file = "./data/temp/Mast.rds")
Endocrine <- readRDS(file = "./data/temp/Endocrine.rds")
Plasma <- readRDS(file = "./data/temp/Plasma.rds")
Schwann <- readRDS(file = "./data/temp/Schwann.rds")

###T_sub_ALL_DimPlot####
pdf("./data/output/T_sub_ALL_DimPlot.pdf",width = 13,height = 6)
DimPlot(Fibroblast, reduction = "umap", split.by = "tech")
DimPlot(B, reduction = "umap", split.by = "tech")
DimPlot(T, reduction = "umap", split.by = "tech")
DimPlot(Ductal, reduction = "umap", split.by = "tech")
DimPlot(Macrophage, reduction = "umap", split.by = "tech")
DimPlot(Endothelial, reduction = "umap", split.by = "tech")
DimPlot(NK, reduction = "umap", split.by = "tech")
DimPlot(Acinar, reduction = "umap", split.by = "tech")
DimPlot(Neutrophil, reduction = "umap", split.by = "tech")
DimPlot(Stellate, reduction = "umap", split.by = "tech")
DimPlot(Mast, reduction = "umap", split.by = "tech")
DimPlot(Endocrine, reduction = "umap", split.by = "tech")
DimPlot(Plasma, reduction = "umap", split.by = "tech")
DimPlot(Schwann, reduction = "umap", split.by = "tech")
dev.off()
#亚群再聚类####
##1. T亚群聚类####
time_start <- Sys.time()
T<-readRDS("./data/temp/T.rds")
#可以不跑 JackStraw，ScoreJackStraw
pbmc <- JackStraw(T, num.replicate = 100,dims = 40)
pbmc <- ScoreJackStraw(pbmc, dims = 1:30)
###T_JackStrawPlot_ElbowPlot####
pdf("./data/output/T_JackStrawPlot_ElbowPlot.pdf",width = 8,height = 6)
JackStrawPlot(pbmc, dims = 1:30)
ElbowPlot(pbmc,ndims = 30)
dev.off()
# find how many 15cluster
hms_neighbor<- FindNeighbors(pbmc, dims = 1:30)
obj <- FindClusters(hms_neighbor, resolution = seq(0.5,1.2,by=0.1))
#resolution設置在0.4-1.2之間,越大clusters越多,查看拐點
###T_clustree####
pdf("./data/output/T_clustree.pdf",width = 8,height = 6)
clustree(obj)
dev.off()
hms_cluster <- FindClusters( hms_neighbor, resolution = 1.2)
head(Idents(hms_cluster), 5)
hms_cluster<- RunUMAP(hms_cluster, dims = 1:30)
###T_DimPlot####
pdf("./data/output/T_sub_DimPlot.pdf",width = 7,height = 6)
DimPlot(hms_cluster, reduction = "umap")
dev.off()
saveRDS(hms_cluster, file = "./data/temp/T_test_1.2.rds")

scRNA.markers <- FindAllMarkers(hms_cluster, 
                                only.pos = TRUE,  #特异性高表达marker
                                min.pct = 0.05, 
                                logfc.threshold = 0.05)
write.table(scRNA.markers,file="./data/temp/TcellMarkers.txt",sep="\t",row.names=F,quote=F)

##挑选每个细胞亚群中特意高表达的20个基因
top20 <- scRNA.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC) 
write.csv(file="./data/temp/top20_cell_markers.csv",top20)
#整理成表格，只显示基因名字
top20_table=unstack(top20, gene ~ cluster)
names(top20_table)=gsub("X","cluster",names(top20_table))
write.csv(file="./data/temp/Ttop20_marker_genes_1.2.csv",top20_table,row.names=F)

pdf("./data/output/T_sub_DimPlot.pdf",width = 7,height = 6)
DimPlot(hms_cluster, reduction = "umap")
dev.off()
#耗时
time_end <- Sys.time()
time_cost <- time_end-time_start
time_cost
# Time difference of 20.48172 mins


##2. T亚群注释####
hms_cluster<-readRDS("./data/temp/T_test_1.2.rds")
DimPlot(hms_cluster, reduction = "umap")

fp1 <- FeaturePlot(hms_cluster,reduction ="umap",features =c("GZMK"))
fp1
fp2 <- VlnPlot(hms_cluster, features = c("MS4A1", "CD79A"))
fp2

#细胞及细胞中基因与RNA数量
slotNames(hms_cluster)
#assay
hms_cluster@assays
dim(hms_cluster@meta.data)
View(hms_cluster@meta.data)

new.cluster.ids <- c("CD4_Tn",   "CD8_Trm",  "CD8_Te",                       "CD4_Tem",                   "CD4_Tem", 
                     "CD4_Treg", "CD8_CTL",  "Epithelial",                   "Macrophages",               "CD4_Treg", 
                     "NKT",     "CD4_Treg", "Macrophages Tumor-associated", "Macrophages Proliferating", "Dendritic cell", 
                     "Dendritic cell", "Acinar cell", "CD4_Th")
names(new.cluster.ids) <- levels(hms_cluster)
hms_cluster_id<- RenameIdents(hms_cluster, new.cluster.ids)
pdf("./data/output/T_hms_cluster_id_DimPlot1.pdf",width =23 ,height = 10)
#未标注有标签 
p1 <- DimPlot(hms_cluster, reduction = "umap", label = TRUE, pt.size = 0.5)
#已标注有标签
p2 <- DimPlot(hms_cluster_id, reduction = "umap", label = TRUE, pt.size = 0.5) 
#无标注有标签
p3 <- DimPlot(hms_cluster_id, reduction = "umap", label = TRUE, pt.size = 0.5)+
  NoLegend()
#已标注无标签
#p4 <- DimPlot(hms_cluster_id, reduction = "umap", label = FALSE, pt.size = 0.5)
g1_treat <- WhichCells(hms_cluster_id, idents = c("CD4_Tn",   "CD4_Th", "CD4_Tem","CD4_Treg", 
                                                  "NKT", "CD8_Trm",  "CD8_CTL",  "CD8_Te"))
# colors <- c("#942d8d","#3cb346","#00abf0","#d75427",
#             "#2e409a","#FB8072","#eeb401")#,"#7BAFDE"
# DimPlot(hms_cluster_id,label=T,cells = g1_treat, cols =colors )
p4 <- DimPlot(hms_cluster_id,label=T,cells = g1_treat)
p1+p2
p3+p4
dev.off()

saveRDS(hms_cluster_id, file = "./data/temp/T_cluster_id_test.rds")
T_cluster_id<-readRDS("./data/temp/T_cluster_id_test.rds")

##3. T亚群注释验证####
hms_cluster<-hms_cluster_id
hms_cluster<-readRDS("./data/temp/T_cluster_id_test.rds")

markers1 <- c(
  "CCR7",  "LTB",# CD4_Tn
  "CCL4", "CCL5", # CD8_Trm 
  "GZMA", "GZMH", # CD8_Te "GZMK", "CCL5",
  "IL7R", "KLRB1",  "AQP3", # CD4_Tem"IL7R",  
  "IL2RA",  "CTLA4", # CD4_Treg  "TIGIT", "BATF","TNFRSF4","TNFRSF18",
  "GZMK","CRTAM",  # CD8_CTL"GZMK", 
  "GZMB","GNLY",# NKT
  "TNFRSF18", "TNFRSF4")   # CD4_Th
dp1 <-DotPlot(subset(hms_cluster,idents=c("CD4_Tn","CD4_Tem","CD8_Te","NKT", "CD4_Th","CD4_Treg","CD8_Trm","CD8_CTL")), #
              features = markers1, 
              dot.scale = 5) + 
  RotatedAxis()

dp1

markers2=c(
  "KRT19",                 # Epithelial
  "GNLY", "GZMB", "KLRD1", "PRF1",# NK cells
  "IL1B", "CXCL8",         # Macrophages Tumor-associated
  "TOP2A","STMN1","MKI67", # Macrophages Proliferating
  "HLA-DRB1", "HLA-DQA2", "HLA-DQA1",  # Dendritic cell
  "PRSS1","CLPS","PRSS2") # Acinar cell"REG1B",
dp2 <-DotPlot(subset(hms_cluster,idents=c("Epithelial", "Macrophages", "Macrophages Tumor-associated", "Macrophages Proliferating", "Dendritic cell", "Dendritic cell", "Acinar cell")), 
              features = markers2, 
              dot.scale = 5) + 
  RotatedAxis()
dp1+dp2
####DotPlot_T亚群聚类_标注验证####
pdf("./data/output/DotPlot_T亚群聚类_标注验证.pdf",width = 13,height = 6)
dp1
dev.off()
###⚠️T亚群注释的并不清晰####
pdf("./data/output/DotPlot_T亚群聚类_标注验证1.pdf",width = 17,height = 6)
dp1+dp2
dev.off()

fp1 <- FeaturePlot(hms_cluster,reduction ="umap",features =c("HBB"))
fp1
fp2 <- VlnPlot(hms_cluster, features = c("MS4A1", "CD79A"))
fp2



## 4. T亚群分组绘图####
###Only T cells####

# Only_T<-subset(hms_cluster_id, idents=c("Naive_CD4_T", "Natural_Killer_T","T_Helper",
#                                         "Effector_Memory_CD4_T","Resident_Memory_CD8_T",
#                                         "Terminally_Exhausted_CD8_T","Regulatory_T",
#                                         "Cytotoxic_CD8_T","Effector_Memory_CD8_T",
#                                         "Pre_Exhausted_CD8_T"))
Only_T <- subset(hms_cluster_id, idents=c("CD4_Tn",   "CD4_Th", "CD4_Tem","CD4_Treg", 
                                          "NKT", "CD8_Trm",  "CD8_CTL",  "CD8_Te"))
DimPlot(Only_T, reduction = "umap", label = TRUE, pt.size = 0.5) 
DimPlot(Only_T, reduction = "umap", label = FALSE, pt.size = 0.5) 
saveRDS(Only_T, file = "./data/temp/Only_T_cluster_id_test.rds")


#⚠️⚠️⚠️这里绘图有瑕疵的这里，绘图gene还要重新选的，出图的基因太少了。
###CD4_T绘图####
# "CD4_Tn", "CD4_Th", "CD4_Tem","CD4_Treg", 

Only_T <- readRDS("./data/temp/Only_T_cluster_id_test.rds")
CD4_T<-subset(Only_T, idents=c("CD4_Tn", "CD4_Th", "CD4_Tem","CD4_Treg"))
saveRDS(CD4_T, file = "./data/temp/CD4_T.rds")

pdf("./data/output/CD4_T_genes11_DimPlot_RidgePlot_DoHeatmap_Heatmap_genes_celltype.pdf",width = 8,height = 6)
DimPlot(CD4_T, reduction = "umap")
RidgePlot(CD4_T, features = c("IL7R", "LTB", "CD40LG", "AQP3", "MAL", "CCR7", "TIMP1", "GPR183", "ANXA1", "CXCR4", "CD4", "CXCL13", "KLRB1", "CSTB"),
          cols = c("green3","orangered"), group.by="tech", ncol = 4) + theme(axis.title.y = element_blank())

markers.to.plot<-c("IL7R", "LTB", "CD40LG", "AQP3", "MAL", "CCR7", "TIMP1", "GPR183", "ANXA1", "CXCR4", "CD4", "CXCL13", "KLRB1", "CSTB")
DoHeatmap(subset(CD4_T,downsample=50000),features = markers.to.plot,size=5)

#expression level in  CD4_T sub cell 
CD4_T_genes11_heatmap<-DotPlot(CD4_T,features = c("IL7R", "LTB", "CD40LG", "AQP3", "MAL", "CCR7", "TIMP1", "GPR183", "ANXA1", "CXCR4", "CD4", "CXCL13", "KLRB1", "CSTB"))+RotatedAxis()
CD4_T_genes11_heatmap
CD4_T_genes11_heatmap<-CD4_T_genes11_heatmap$data
write.csv(CD4_T_genes11_heatmap,"./data/temp/CD4_T_genes11_heatmap.csv")

#expression level in each patients
CD4_T_celltype_heatmap<-DoHeatmap(subset(CD4_T,downsample=50000),features = markers.to.plot,size=5,group.by="celltype")
CD4_T_celltype_heatmap
CD4_T_celltype_heatmap<-CD4_T_celltype_heatmap$data
write.table(CD4_T_celltype_heatmap,"./data/temp/CD4_T_celltype_heatmap.csv")
dev.off()

###CD8_T绘图####

# "CD4_Tn",   "CD4_Th", "CD4_Tem","CD4_Treg", 
CD8_T<-subset(Only_T, idents=c("CD8_Trm",  "CD8_CTL",  "CD8_Te"))
saveRDS(CD8_T, file = "./data/temp/CD8_T.rds")

pdf("./data/output/CD8_T_genes11_DimPlot_RidgePlot_DoHeatmap_Heatmap_genes_celltype.pdf",width = 8,height = 6)
DimPlot(CD8_T, reduction = "umap")
RidgePlot(CD8_T, features = c("CD8A","CD8B","CD38","CD69","ENTPD1","GZMA","GZMH","MYO1F","SYNE1","TSC22D3","XCL2"),
          cols = c("green3","orangered"), group.by="tech", ncol = 4) + theme(axis.title.y = element_blank())

markers.to.plot<-c("CD8A","CD8B","CD38","CD69","ENTPD1","GZMA","GZMH","MYO1F","SYNE1","TSC22D3","XCL2")
DoHeatmap(subset(CD8_T,downsample=50000),features = markers.to.plot,size=5)

#expression level in  CD8_T sub cell 
CD8_T_genes11_heatmap<-DotPlot(CD8_T,features = c("CD8A","CD8B","CD38","CD69","ENTPD1",
                                                     "GZMA","GZMH","MYO1F","SYNE1","TSC22D3","XCL2"))+RotatedAxis()
CD8_T_genes11_heatmap
CD8_T_genes11_heatmap<-CD8_T_genes11_heatmap$data
write.csv(genes11_heatmap,"./data/temp/CD8_T_genes11_heatmap.csv")

#expression level in each patients
CD8_T_celltype_heatmap<-DoHeatmap(subset(CD8_T,downsample=50000),features = markers.to.plot,size=5,group.by="celltype")
CD8_T_celltype_heatmap
CD8_T_celltype_heatmap<-celltype_heatmap$data
write.table(celltype_heatmap,"./data/temp/CD8_T_celltype_heatmap.csv")
dev.off()


###Only_T绘图####
Only_T <- readRDS(file="./data/temp/Only_T_cluster_id_test.rds") 

# 修改ident的顺序，解决后期作图排序
levels(x = Only_T)
levels(x = Only_T) <- c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", "NKT", "CD8_CTL", "CD8_Trm", "CD8_Te")
levels(x = Only_T)


markerdata$subtype <- factor(x=markerdata$celltype,
                              levels = c("Endothelial","Fibroblast","Epithelial","Immune","Other"))



p1 <- DimPlot(Only_T, reduction = "umap", label = TRUE, pt.size = 0.5) 
p1
p2 <- RidgePlot(Only_T, features = #c("CD8A","CD8B","CD38","CD69","ENTPD1", "GZMA","GZMH","MYO1F","SYNE1","TSC22D3","XCL2"))+RotatedAxis()
                                    c("CCR7",  "LTB",# CD4_Tn
                                      "TNFRSF18", "TNFRSF4"  , # CD4_Th
                                      "IL2RA",  "CTLA4", # CD4_Treg  "TIGIT", "BATF","TNFRSF4","TNFRSF18",
                                      "IL7R", "KLRB1",  "AQP3", # CD4_Tem"IL7R", 
                                      "GZMB","GNLY",# NKT
                                      "GZMK","CRTAM",  # CD8_CTL"GZMK", 
                                      "CCL4", "CCL5", # CD8_Trm 
                                      "GZMA", "GZMH"), # CD8_Te "GZMK", "CCL5", 
            cols = c("green3","orangered"), group.by="tech", ncol = 4) + theme(axis.title.y = element_blank())
p2

#T细胞热图
markers.to.plot <- #c("CD8A","CD8B","CD38","CD69","ENTPD1", "GZMA","GZMH","MYO1F","SYNE1","TSC22D3","XCL2"))+RotatedAxis()
                    c("CCR7",  "LTB",# CD4_Tn
                      "TNFRSF18", "TNFRSF4"  , # CD4_Th
                      "IL2RA",  "CTLA4", # CD4_Treg  "TIGIT", "BATF","TNFRSF4","TNFRSF18",
                      "IL7R", "KLRB1",  "AQP3", # CD4_Tem"IL7R", 
                      "GZMB","GNLY",# NKT
                      "GZMK","CRTAM",  # CD8_CTL"GZMK", 
                      "CCL4", "CCL5", # CD8_Trm 
                      "GZMH", "GZMA") # CD8_Te "GZMK", "CCL5", 
p3 <- DoHeatmap(subset(Only_T,downsample=50000),
                features = markers.to.plot,
                size=5,
                group.bar.height=0.03)+
  scale_fill_gradientn(colors = c("white","grey","firebrick3"))
p3

Only_T_genes11_heatmap<-DotPlot(Only_T,features = #c("CD8A","CD8B","CD38","CD69","ENTPD1", "GZMA","GZMH","MYO1F","SYNE1","TSC22D3","XCL2"))+RotatedAxis()
                                                  c("CCR7",  "LTB",# CD4_Tn
                                                    "TNFRSF18", "TNFRSF4"  , # CD4_Th
                                                    "IL2RA",  "CTLA4", # CD4_Treg  "TIGIT", "BATF","TNFRSF4","TNFRSF18",
                                                    "IL7R", "KLRB1",  "AQP3", # CD4_Tem"IL7R", 
                                                    "GZMB","GNLY",# NKT
                                                    "GZMK","CRTAM",  # CD8_CTL"GZMK", 
                                                    "CCL4", "CCL5", # CD8_Trm 
                                                    "GZMA", "GZMH")) # CD8_Te "GZMK", "CCL5", 
Only_T_genes11_heatmap
Only_T_genes11_heatmap<-Only_T_genes11_heatmap$data
write.table(Only_T_genes11_heatmap,"./data/temp/Only_T_genes11_heatmap.csv")

markers.to.plot<-#c("CD8A","CD8B","CD38","CD69","ENTPD1", "GZMA","GZMH","MYO1F","SYNE1","TSC22D3","XCL2"))+RotatedAxis()
                  c("CCR7",  "LTB",# CD4_Tn
                    "TNFRSF18", "TNFRSF4"  , # CD4_Th
                    "IL2RA",  "CTLA4", # CD4_Treg  "TIGIT", "BATF","TNFRSF4","TNFRSF18",
                    "IL7R", "KLRB1",  "AQP3", # CD4_Tem"IL7R", 
                    "GZMB","GNLY",# NKT
                    "GZMK","CRTAM",  # CD8_CTL"GZMK", 
                    "CCL4", "CCL5", # CD8_Trm 
                    "GZMA", "GZMH") # CD8_Te "GZMK", "CCL5", 
p4<-DoHeatmap(subset(Only_T,downsample=50000),
              features = markers.to.plot,
              size=5,
              group.by="celltype")+
  scale_fill_gradientn(colors = c("white","grey","firebrick3"))
p4
Only_T_celltype_heatmap<-p4$data
write.table(Only_T_celltype_heatmap,"./data/temp/Only_T_celltype_heatmap.csv")

pdf("./data/output/Only_T_genes11_DimPlot_RidgePlot_DoHeatmap_Heatmap_genes_celltype.pdf",width = 16,height = 13)
p1
p2
p3
p4
dev.off()











